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A method for estimating the axis of reflectional symmetry of an 
image f{x,y) on the unit disc D = {(x,y) -x^ + <1} is proposed, 
given that noisy data of f{x,y) are observed on a discrete grid of 
edge width A. Our estimation procedure is based on minimizing over 
P £ [0,7r) the L2 distance between empirical versions of / and rpf, 
the image of / after reflection at the axis along (cos /3, sin Here, 
/ and Tjif axe estimated using truncated radial series of the Zernike 
type. The inherent symmetry properties of the Zernike functions re- 
sult in a particularly simple estimation procedure for /3. It is shown 
that the estimate /3 converges at the parametric rate for images 
/ of bounded variation. Further, we establish asymptotic normality 
of /3 if / is Lipschitz continuous. The method is applied to calibrat- 
ing the point spread function (PSF) for the deconvolution of images 
from confocal microscopy. For various reasons the PSF characterizing 
the problem may not be rotationally invariant but rather only reflec- 
tion symmetric with respect to two orthogonal axes. For an image 
of a bead acquired by a confocal laser scanning microscope (Leica 
TCS), these axes are estimated and corresponding confidence inter- 
vals are constructed. They turn out to be close to the coordinate axes 
of the imaging device. As cause for deviation from rotational invari- 
ance, this indicates some slight misalignment of the optical system 
or anisotropy of the immersion medium rather than some irregular 
shape of the bead. In an extensive simulation study, we show that us- 
ing a symmetrized version of the observed PSF significantly improves 
the subsequent reconstruction process of the target image. 
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1. Introduction. The fundamental concept of symmetry of physical and 
biological objects has been thoroughly studied for a long time; cf., e.g., Con- 
way, Burgiel and Goodman-Strauss (2008). In particular, symmetry plays 
an important role in image analysis and understanding and finds direct ap- 
plications in object recognition, robotics, image animation and image com- 
pression; see Liu, Collins and Tsin (2004) for an overview of the subject 
of symmetry and related issues. The problem of detecting and measuring 
object symmetries has been tackled in the image processing and pattern 
analysis literature since the original works of Atallah (1985) and Friedberg 
(1986). For a comprehensive review of the literature see Liu, Collins and Tsin 
(2004) and Bissantz, Holzmann and Pawlak (2009). The role of symmetry 
in statistical inference is discussed in Viana (2008). 

In this paper we propose an estimation procedure for the angle /3 of the di- 
rection (cos /3, sin /3) of the axis of reflectional symmetry of an image function 
/ from which discrete, noisy observations are available. The observations are 
taken on a grid of edge width A, and the noise is modeled by stochastic er- 
rors. Existing methods either do not allow any noise or treat the effect of 
noise only empirically by simulations. Thus, to the best of our knowledge, 
our approach is the first which treats refiection symmetry estimation from 
a statistical point of view as a semiparametric estimation problem. Specifi- 
cally, we show that for image functions / of bounded variation the estimate 
13 converges at a rate of and, further, for Lipschitz continuous / we have 
asymptotic normality, which allows us to construct asymptotic confidence 
intervals for /?. 

The estimation procedure is based on minimizing over /3 € [0,7r) the L2 
distance between empirical versions of / and r^/, the image of / after re- 
flection at the axis along (cos/3, sin/3). Here, / and r^/ are estimated using 
truncated Zernike function expansions. The inherent symmetry properties 
of the Zernike functions yield a particularly simple estimation procedure for 
p. In the recent related papers [cf. Kim and Kim (1999) and Revaud, Lavoue 
and Baskurt (2008)], methods for estimating the rotation angle of an im- 
age invariant under a certain rotation, which also make use of the Zernike 
moments, have been proposed. However, the authors do not study any con- 
vergence aspects of the algorithms and confine their discussion to noise-free 
images. 

Our methodology is applied to calibrating the point spread function (PSF) 
of a microscope in confocal microscopy. The PSF describes the blurring effect 
of the imaging process. Typical smoothing scales are of order wlOO nm, 
which often is of similar order as the size of relevant structures in the target 
object. Hence, an exact knowledge of the PSF is essential to properly adjust 
(i.e., deconvolve) the observed image to recover the image of the target 
object. 
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A theoretical PSF may be computed from the optical properties of the 
microscope, it is rotationally invariant for a rotationally symmetric optical 
system. However, the true (empirical) PSF can deviate substantially from 
its theoretical shape, and is no longer rotationally invariant. Therefore, the 
PSF is estimated from images of point-like objects with known form. Since 
this process involves rather dim images, it is worthwhile to use additional 
information on the PSF to improve on its reconstruction. 

Often, the empirical PSF is still expected to be reflection symmetric with 
respect to two (unknown) orthogonal axes, for example, if the detector plane 
is not in perfect agreement with the focal plane of the microscope; cf. Lehr, 
Sibarita and Chassery (1998) and Pankajakshan et al. (2008). Therefore, for 
an image of a bead acquired by a confocal laser scanning microscope (Leica 
TCS), in Bissantz, Holzmann and Pawlak (2009) we used hypotheses tests 
to assess rotational invariance as well as invariance under a rotation by vr 
(which is a consequence of invariance under reflections by two orthogonal 
axes) for the empirical PSF. While (for bead 2) rotational invariance was 
rejected, invariance under a rotation by vr (and hence reflection symmetry) 
was not rejected at the level of 5%. 

Here, we estimate the axes of reflectional symmetry of the PSF and con- 
struct the corresponding confidence intervals. It turns out that the axes are 
very close to the coordinate axes of the imaging device. This indicates that 
the reason for the PSF to deviate from rotational invariance appears to be 
some (slight) misalignment of the optical system or anisotropy of the immer- 
sion medium used for object preparation rather than some random deviation 
from sphericity of the bead used to image the PSF. 

Further, we propose to reduce the noise level in the PSF by a factor of 
2 by averaging along the estimated axes. To investigate the practical merit 
of this strategy for recovery of a target image, we use a two-step simula- 
tion study. First, the PSF is estimated by four different methods, then the 
estimated PSFs are used for subsequent recovery of the target image, and 
the accuracy of these reconstructions are compared. For the PSF we use a 
simple nonparametric estimate of the PSF as well as a symmetrized ver- 
sion, together with correctly specified and slightly misspecified parametric 
models. It turns out that while the correctly specified parametric model per- 
forms best for recovering the target image, symmetrizing the nonparametric 
estimate greatly improves its performance, even beyond that of the slightly 
misspecified parametric model. 

The paper is organized as follows. In Section 2 we introduce the theoretical 
Zernike moments and give their basic invariance properties. Further, we dis- 
cuss how to estimate the moments from data generated by our observational 
model. In Section 3 we propose the estimation procedure for the angle /3 of 
the direction (cos/?, sin /3) of the axis of reflectional symmetry of the image 
function /, and discuss its statistical properties. This includes the issue of 
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uniqueness as well as consistency, rate of convergence and asymptotic distri- 
bution of the estimate. Section 4 contains simulation studies concerning the 
finite sample properties of the estimator. In Section 5 we discuss reflection 
symmetry properties of an observed PSF of a confocal laser scanning mi- 
croscope (Leica TCS). Further, in a simulation we show how incorporating 
reflection symmetry into a simple nonparametric estimate of the PSF signif- 
icantly improves its properties in the image reconstruction process. Section 
6 gives some concluding remarks, while technical proofs can be found in the 
supplementary material in Bissantz, Holzmann and Pawlak (2010). 

2. The Zernike orthogonal basis and image reconstruction. Zernike func- 
tions, introduced as an orthogonal and rotationally invariant basis of poly- 
nomials on the disc in Zernike (1934), and their corresponding moments 
have been used extensively in image analysis and pattern recognition; see 
Bailey and Srinath (1996), Khotanzad and Hong (1990) and Mukundan and 
Ramakrishnan (1998). The Zernike basis has also been employed as an im- 
portant tool for the statistical inference concerning the inverse problem of 
positron emission tomography [cf. Jones and Silverman (1989) and John- 
stone and Silverman (1990)] and PSF estimation in fluorescence microscopy 
[cf. Dieter len et al. (2004); Dieter len et al. (2008)]. 

2.1. Zernike polynomials. In the following we identify two-dimensional 
space M? with the complex plane C via {x,y) ^ x + iy, where i is the imagi- 
nary unit. In particular, e*^ is the unit vector (cos /3, sin /3) at angle /? to the 
X axis. 

Now, the (complex) Zernike orthogonal polynomials are given by Vpq{x, y) = 
Rpq{p)e'^'^^ , {x,y) G D, where p = yx^ + y^, 9 = arctan(y/x) and Rpq{p) is 
the radial Zernike polynomial given explicitly by 

... ''1^" {-mp-iw~-^ 

h ^K(p + kl)/2-0!((p-kl)/2-0!' 

The indices (p, q) have to satisfy p > 0, \q\ < p, and p — \q\ has to be even. 
We will call such pairs (p, q) admissible. The Zernike polynomials satisfy the 
following orthogonality relation over the unit disc D: 

Vpqix, y)V*,q, (X, y) dx dy = Tl/{p + l)5pp,5qq>, 

Id 

where * denotes complex conjugation and 5ppi is the Kronecker delta. This 
implies that 

(1) ll^pgf = V(P + 1)=^P, 

where || • || is the norm on L2{D). In Bhatia and Wolf (1954), the Zernike 
polynomials are characterized by a certain uniqueness property, among oth- 
ers, invariant polynomials defined on D. 
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2.2. Function approximation. Since the family {Vpq{x,y)} for admissible 
{p,q) forms a complete and orthogonal system in L2{D), we can expand a 
function / € L2{D) into a series of the Zernike polynomials, that is. 



oo p 



(2) fix,y) = Y^ Yl n;'Ap,{f)Vp,{x,y), 

p=0 q=—p 

where here and throughout the paper the summation is taken over admissible 
pairs {p,q). Thus, the Fourier coefficients {Apg{f)} (often referred to as the 
Zernike moments) uniquely characterize the image function /. The norming 
factor arises due to (1), and the Zernike moment Apq{f) is defined by 

Apqi f) = jj^fix,y)Vpg{x,y)dx dy. 
Owing to Parseval's formula, we have that for / G L2{D) 

oo p 

(3) ll/f = EE 

p=0 q=—p 

Let us introduce the notation f{p,9) = f{pcos9,psm9) for a function / € 
L2{D). Then by using polar coordinates we obtain 



(4) 



\qU) = 27r / Cg{p,f)Rpg{p)pdp, 
Jo 

cq{pj) = ^fj f{p,e)e-"iUe. 



2.3. Image reconstruction. We assume that the data are observed on a 
symmetric square grid of edge width A, that is, Xi — Xi-i = yi — yi-i = A 
and Xi = —Xm-i+i, yi = —ym-i+i, so that {xi,yj) is the center of the pixel 
Ilij = [xi — ^ ,Xi + ^] X [yj — ^,yj + Note that m corresponds to 2/ A. 
For / G L2{D) we shall assume the following observational model: 

(5) Zij = f{xi,yj) + eij, {xi,yj) e D,l <i,j <m, 

where the noise process {eij} is an i.i.d. random sequence with zero mean, 
finite variance Eefj = cj^ and finite fourth moment, so that Zij is the datum 
associated with pixel Hij. Note that along the boundary of the disc, some 
lattice squares are included (if their center is in D) and some are excluded. 
When reconstructing /, this gives rise to an additional error, called geometric 
error in Pawlak and Liao (2002). This error can be quantified by using the 
celebrated problem in analytic number theory referred to as lattice points 
of the circle. In applications, the datum Zjj might also correspond to the 
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average of / over the pixel Iljj rather than its value at the center, in such 
cases we assume negligible variation of / over Uij. 

In the following we need to work with a discretized version of the Zernike 
moments. Consider weights Wpq{xi,yj) of the form 

(6) Wpq{xi,yj) = jj^ V*g{x,y)dxdy or Wpg{xi,yj) = A'^V*g{xi,yj). 
Using either version in (6), we estimate the Zernike moment Apq(f) by 

(7) Apq= ^ Wpq{Xi,yj)Zij. 

{xi,yj)&D 

For efficient methods for computing the Zernike moments Apq^ see, for ex- 
ample, Amayeh et al. (2005). Instead of the uniform weights, one could also 
use a more sophisticated quadrature rule, particularly if sharp features of / 
are expected. 

3. Reflection estimation. First we investigate the effect that reflecting an 
image function / has on its Zernike moments. Suppose that / is reflected at 
a line along the direction e*^, /3 G [0,vr), and denote the reflected function by 

Tpf . Then one easily shows that {Tpf){p, 6) = /(p, 2/3 — 9) and, consequently, 
using (4), 

(8) Apq{Tpf) = e-^^^^Ap,.q{f). 

Consider the following assumption. 

Assumption 1. Suppose that / G L2{D) is invariant under some unique 
reflection r^. . 

Indeed, the composition of two reflections along lines e*°^ and e*°^ is a 
rotation with angle 2 (02 — cti). Thus, / is invariant under a unique reflection 
if and only if / is invariant under some reflection and if / is not invariant 
under any rotation. 

3.1. Contrast functions. Our method for estimating f3* is based on the 
expansion (3) and the invariance property of Zernike moments expressed by 
the formula in (8). We set 

00 p 

(9) M(/3,/) = ||/-Wf = E^p' E \Apq{f)-e-'^^^Ap,.q{f)\\ 

p=o q=—p 
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Evidently, under Assumption 1 the angle /3* is the unique zero of the func- 
tion M{P,f). Writing Apg{f) = \Apg{f)\e''-'''-f^ and noting that Apgif) = 
Ap^^q{f)* , we calculate 

oo p 

(10) M(/3, /) = ^n;i ^4|^p,(/)|2(l - cos(2rp,(/) + 2q(3)), 

where the sums are taken over admissible pairs {p,q). Therefore, /3* is also 
uniquely characterized by the condition cos(2rpq(/) + 2q/3*) = 1 or by re- 
quiring 

(11) rp,{f)eqP*+7rZ 

for all p, q with Apg(f) / 0. 

Thus, a natural way to estimate /3* is to first estimate a truncated version 
of the series defining M(/3,/), and then define an estimate of /3* as the 
minimizer of this estimated contrast function. We first show that suitably 
truncated versions of M(/3,/) still uniquely determine /3*. For a fixed set 

TV p 

MM{PJ)=Y.n-^ \Ap,{f) - e-'^'^P Ap,.,{f)\\ 
p=0 q=—p 

This is the truncated counterpart of the series in (9). Evidently, Mi^(f3* , f) = 
for all N under Assumption 1. We shah call M(/3, /) and Mn{P, f) as well 
as their empirical version below contrast functions. 

Theorem 1. Suppose that f satisfies Assumption 1. Then for suffi- 
ciently large N = N(f), /3* is the unique zero of the truncated contrast func- 
tions M]\f{f3, f). 

The proof of Theorem 1, given in the supplementary material in Bissantz, 
Holzmann and Pawlak (2010), reveals that in order to uniquely determine 
the direction of the reflection axis e^^ as the zero of the function M^{f3, f) 
one has to choose A^ so large such that the sum defining Mi\f{/3,f) contains 
nonzero Apg^s for which the greatest common divisor (gcd) of the g's is 
1. Thus, we can choose A^ as the smallest value such that j4pjqj(/) 7^0, 
•■ - ^-^p^^rl/) / O, for Pi < A^,i = l,...,r, with gcd(gi, . . . , g^) = 1- 

In practice, Mat (/?,/) and hence an appropriate value for A^ still has to 
be estimated. One could test sufficiently many moments to be nonzero, how- 
ever, we prefer to choose A^ for appropriate estimation of / in the resulting 
truncated Zernike series estimate; see Section 5. Apart from its theoretical 
value. Theorem 1 implies that even in dim images occurring, for example, 
in PSF estimation in Section 5, where only few Zernike moments may be 
properly estimated, it is still possible to identify and estimate the symmetry 
axis. 
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3.2. Estimation. For estimation purposes, we first estimate tlie contrast 
functions M]\f{f3,f) by 

N p 

p=o g=-p 

where we write Apg = {Apgle^^pi . Then we define the estimator of /3* as 

/9A,Ar = arg min MmW)- 

/3G[0,7r) 

The estimate /3a, Af depends on the grid size A and, more importantly, on 
the truncation parameter N . Note that although Mj\f{l3*) = 0, M7v(/3a,7v) 
will be positive a.s. due to noise. 

Remark 1 . The estimated contrast function Mjv (/?) is simply the squared 
L2 distance between the Zernike estimate given in polar coordinates by 

_ N 

ip,tl) 

and its reflected version r^/. Note that this is achieved by a special property 
of the Zernike polynomials, namely, the set of Zernike polynomials used in 
the estimate / remains invariant under reflection. As suggested by a referee, 
an estimate similar to f3A,N would be obtained by estimating the coefficients 
Ap^g in 

N 

(12) f{p, 0) - T^fip, e) = J2 V(^P^i - Ap,-,e-"''')Vp,M 0) 

by least squares for each fixed /3, and then choosing the (3 with minimal 
RSS. While our approach is somewhat simpler since we estimate the Ap^q 
before imposing symmetry (and thus independently of /3), this approach 
could potentially be placed into a likelihood or Bayesian framework as in 
Pankajakshan et al. (2008). 

The next result states uniform convergence in probability of the estimated 
contrast function M]\f{l3) to Mj^{f3). This is also used in order to obtain the 
consistency of /3a, at for f3* . 

Theorem 2. For each fixed N, as A — > 0, 

(13) sup |Mjv(/3)-/V/;v(/3,/)|^0 (P), 

/3e[0,7r) 

where (P) denotes convergence in probability. 
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Note that in Theorem 2, / need not be reflection invariant. The next 
theorem gives the consistency of Pa,n as A — )• as weh as its parametric 
A-rate of convergence. In ah the results that follow we choose the truncation 
parameter N according to the prescription established in Theorem 1, that 
is, we require that should be selected in such a way that /3* is the unique 
minimizer of Mi^{(3,f). We will refer to such a value of N as "sufficiently 
large." 

Theorem 3. Suppose that f G L2{D) is a function of bounded variation 
and satisfies Assumption 1. Then for sufficiently large (hut fixed) N = N{f), 
we have that, as A — )■ 0, 

(14) |/3A,iV-/3*|=Op(A). 

Next we establish asymptotic normality for the estimate (3a,n- In order 
for the bias term of the estimated Zernike coefficient to be negligible, we 
require that the image function / is Lipschitz continuous. 

Theorem 4. Suppose that f is Lipschitz continuous and satisfies As- 
sumption 1. Then for sufficiently large (hut fixed) N = N{f), we have that, 
as A 0, 

(15) A-.(;3,,„_;3.)4„(o,^j^), 

where 

N p 

(16) M'j(,i(3*J) = 16|^p.(/)Pg'. 

Theorem 4 can be used to construct an asymptotic confidence interval 
for /3*. To this end, we need an estimate of the asymptotic variance in the 
normal limit (15). We may estimate M^(/3*, /) directly by using (16) simply 
by replacing Apq{f) by Apq. However, this may result in underestimation 
of the asymptotic variance, and therefore plugging /3a, Af into the second 
derivative of Miy{/3), 

N p 

M'M = J^^;' 16|ip,|2g2cos(2rp, + 2(^/3), 

p=0 g=0 

should generally be preferred. Call either estimate M^. Further, we need to 
estimate the error variance o"^. To this end, one could use the residuals from 
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the fitted truncated Zernike series. We prefer to use a difference estimate of 
the form 

^^^^ ^' = cTa) ^ \m,-z,^,,f + {z,,-z,,+,f), 

which does not rely on the same underlying regression estimate. Here the 
sum is taken over all {xi,yj) € D where € D and € D, 

and C(A) is the number of terms in this restricted sum. One can show that 
if / is Lipschitz continuous, then = Op (A). For detailed information 

on difference-based estimators in higher dimensions see Munk et al. (2005). 

Using these estimates, we obtain the following confidence interval with 
nominal level a for P*: 

where ui-a is the 1 — o-quantile of the standard normal distribution. 

Remark 2. If in Theorem 4 we only assume that / G L2{D) is a func- 
tion of bounded variation, then the bias is also of order A, and we get an 
asymptotic offset (i.e., a limiting normal law with nonzero mean) in (15). 

Remark 3. If the image / is not reflection invariant, the estimator /3a, 7V 
may still converge to a certain parameter value (3^, which is determined 
by minimizing the L2-distance ||/ — t^/P- Then (/ + t^»/)/2 is the best 
reflection-symmetric approximation (in the L2 sense) to the original image 
/. However, since (3^ is no longer a zero of the contrast function M{(3,f), 
Theorem 1 does not hold, and in order to achieve consistent estimation 
theoretically, one requires that N ^ 00. 

Remark 4. Suppose that / is reflection invariant but is also invariant 
under some discrete rotation group. Then there will be a minimal angle 
a = 2'K /d for some d € N, under rotation of which / is invariant. If we use 
the estimator /3a, Af in such a situation, then a unique reflection axis will be 
between and a, and one should use the minimizer of M]\f{(3) in the interval 
[0,a) rather than in [0,7r). 

4. Finite sample performance. 

4.1. Target functions and the shape of their contrast functions. In this 
section we discuss the results of a simulation study of the proposed estima- 
tion method for the angle /3 of reflectional symmetry. We performed sim- 
ulations with three target functions, which are given in polar coordinates 
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by 



fi{p,0) 

f2{p,e) 
h{p,e) 



ci • X • (1 - /)) • (sin(y + \/ x'^ + y^) + sin(-y + y'x^^+y^)), 

C2 . p . (1 - . (ecos(e)/0.02 ^ gco.(9+0.6)/0.02 

_|_ gCos(6l-0.3+7r)/0.02 _|_ gCos(0+O.9+7r)/O.O2-j 
C3 . p . (1 - p) . (e™«(^)/0-2 + ecos(e+0.9)/0.2 ^ g . gcos((?-1.7)/0.2^^ 



where x = pcos(^), y = psm{9), and ci,C2,C3 are normalization constants 
such that the squared functions all integrate to one on the unit disc. Figures 
1-3 show the target functions without noise and with Gaussian noise, where 
the signal-to-noise ratio, defined as the ratio between the peak values of 
the respective target function /i,/2,/3 and the standard deviation of the 
noise a, is 5. Figure 4 again shows /i but with a signal-to- noise ratio of 
16.7. Note that the functions fi and /2 are reflection symmetric, whereas 
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Fig. 1. Reflection symmetric function fi without noise, with Gaussian noise, 
and AIr{l3) (full curve) and Mii^fi) (dashed curve). Parameters are n = 25 and 
signal-to-noise-ratio = 5. The vertical line indicates the direction of reflection symmetry 
in the true image. 
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Fig. 2. Reflection symmetric function f2 without noise, with Gaussian noise, 
and MrdS) (full curve) and Mt^P) (dashed curve). Parameters are n — 25 and 
signal-to-noise-ratio = 5. The vertical line indicates the direction of reflection symmetry 
in the true image. 



/s is not. Moreover, in all cases we have used regularization parameters N 
chosen according to the selection rule described in Bissantz, Holzmann and 
Pawlak (2009) (a stochastic analogue of the numerical discrepancy principle 
for parameter selection in inverse problems). The fact that /i and /2, in 
contrast to /s, are reflection symmetric is clearly expressed in the shape of 
the associated contrast functions. Indeed, Mj(/3) is far above zero for f^, in 
contrast to the case of /i and /2, where Mj(^/3) reaches a minimum close to 
zero for noisy data. However, we note that even for /a there still exists a well- 
defined minimum of the contrast function Mj{f3). The right panel in Figure 
3 shows a reflection symmetric version of /a, which has been generated by 
adding a version of /s mirrored w.r.t. the axis given by the direction of the 
minimum Mj(P). 

4.2. Simulated distributions of estimated directions (3. In the second part 
we have simulated the distribution of /3, determined as the minimum of 




Fig. 3. The function /g (which is not reflection symmetric) without noise, with Gaussian 
noise, and Mt^P) (full curve), M-j^jS) (dashed curve) and a symmetrized version of f^. 
Parameters are n = 25 and signal-to-noise-ratio = 5. 



M]\f{l3), for a range of values for the parameters n and the signal-to-noise 
ratio s/n. Figures 5 and 6 show density plots of the simulated distributions 
together with normal limits. For the reflection symmetric functions /i and 
/2 we compare the simulated distributions to their asymptotic counterparts 
according to (15). Even for images of moderate size such as the unit circle 
in the square image with edge length (2m + 1) = 51 pixels, the simulated 
distributions are already close to their asymptotic limit. 



5. Calibrating the PSF in confocal microscopy. 



5.1. Assessing reflectional symmetry of the PSF. In this section we use 
the contrast function to estimate the axes of reflection symmetry in an image 
of the point-spread function in confocal fluorescence microscopic imaging. 
Here one observes count data representing observed pixel-integrated image 
intensities on a two-dimensional (or three-dimensional) equidistant grid of 
pixels. We consider the two-dimensional case, where the observations are 
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pi/4 pi/2 3pi/4 pi 

P 



Fig. 4. The reflection symmetric function fi without noise, with Gaussian noise, 
and AIr{l3) (full curve) and Mt{I3) (dashed curve). Parameters are n — 25 and 
signal-to-noise-ratio = 16.7. The vertical line indicates the direction of reflection symmetry 
in the true image. 

Zij = {K-f){xi,yj) + Eij, with 

(19) {K-f){x,y) = k*-f{x,y) = k{x - h.y - t2)^{ti,t2) dtidt2, 

and where "*" represents the convolution of the "true" image 7 G with 
the so-cahed point-spread-function (PSF) k ^ of the microscope. The 
standard model for the distribution of the photon count data Zij is that 
Zi^j is Poisson with the mean {K^){xi,yj), all independent. 

The PSF represents the image of a point-source observed by the micro- 
scope and describes the blurring effect of the imaging process. As discussed 
in the introduction, the PSF is typically estimated by observing a point- 
like object (called bead) of known form. Figure 7 (right) shows the image 
of a bead under a Leica TCS confocal laser scanning microscope. The (ob- 
served) empirical PSF is typically no longer rotationally invariant, but it 
often remains reflection symmetric under two (unknown) orthogonal axes, 
even if, for example, the detector plane was not in perfect agreement with 



n=25,s/n=5 



n=50,s/n=5 



n=100,s/n=5 




qI : ^ : 1 Ql ...^ ^ ^ ^ 1 0^ ' ' ' ' 

0.285 0.29 0.295 0.3 0.305 0.31 0.315 0.294 0.296 0.298 0.3 0.302 0.304 0.306 0.294 0.296 0.298 0.3 0.302 0.304 0.306 

P P P 

Fig. 5. Simulated (solid curve) and asymptotic (dashed curve) distributions of Pa.n for fi. The variance of the asymptotic distributions 
is given as (cf. Theorem 3). The parameter N was chosen as 7,8,12 (first row, left to right) and 8,12,12 (second row, left to 

right). 
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Fig. 6. Simulated (solid curve) and asymptotic (dashed curve) distributions of /3a,n for 
/2 (left panels) and fs (right panels). The variance of the asymptotic distributions is given 
as ^"J^,-^ (cf. Theorem 3). The parameter N was chosen as 7, except for the lower left 
where it was chosen as 8. 



the focal plane of the microscope; cf. Lehr, Sibarita and Chassery (1998) 
and Pankajakshan et al. (2008). 

In Bissantz, Holzmann and Pawlak (2009) we applied tests both for ro- 
tational invariance and for invariance under a rotation by tt (which is an 
immediate consequence of reflection symmetry w.r.t. two orthogonal axes) 
to the observed PSF in the image bead (Figure 7). It turned out that ro- 
tational invariance could be rejected at a 5% level, but invariance under a 
rotation by vr was not rejected. 

For a deeper investigation, we now apply our methodology to estimate 
the (orthogonal) axes of reflection symmetry. The data from fluorescence 
microscopic imaging in general is distributed (approximately) according to 
a Poisson distribution with expectation given by the respective image in- 
tensity. Hence, the noise is not homoscedastic as required by model (5). As 
suggested by a referee, we use the (variance stabilizing) Anscombe transform 
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Fig. 7. Contrast function of image bead (left), image bead with superposed estimated re- 
flection axis and the associated asymptotic confidence interval with nominal level 95% 
(middle), and bead after averaging along two estimated axes of reflectional symmetry 
(right). Bead was acquired during two observation runs of HeLa cervix carcinoma cells 
with a Leica TCS laser scanning fluorescence microscope. Here, we used N = 4. 

[Anscombe (1948)]. Note that reflection symmetry is preserved in this pro- 
cess. Further, fohowing Remark 4, we restrict the range of /3 to [7r/4, 3-7r/4], 
which yields an estimated angle and associated 95% confldence interval of 
/3 = 1.54 lb 0.07. The truncation parameter was selected as = 4 by the 
method described in Bissantz, Holzmann and Pawlak (2009). Using the un- 
transformed data and ignoring heteroscedasticity yields quite similar results 
(/3 = 1.53 lb 0.08), thus, heteroscedasticity appears to be a minor problem in 
this context. 

Figure 7 (left and middle) shows the contrast function (of the untrans- 
formed data with = 4) and the image with superposed estimated reflection 
axis (/3 ?a 1.53). The optical axis appears to be rather close to the coordi- 
nate axis of the image. In particular, the coordinate axis is covered by the 
associated 95%-nominal level confidence interval for f3 [cf. (18)]. This indi- 
cates that the reason for a PSF which is not rotationally invariant appears 
to be some (slight) misalignment of the optical system or anisotropy of the 
immersion medium used for object preparation rather than some random 
deviation from sphericity of the bead used to image the PSF. In Figure 7 
(right), we plot the PSF after averaging along two estimated axes of reflec- 
tional symmetry. 

5.2. Performance of symmetrized PSF estimates for image reconstruction. 
In this section we discuss the results from an extensive simulation study 
in which we investigate the potential beneflt of incorporating symmetry 
information into PSF estimates. 

We shall compare the performance of several models for the PSF for 
subsequent image reconstruction in a two-step simulation procedure which 
mimics the observational process in confocal microscopy. In the first step, 
we generate an image of a point-like object and use it to estimate the PSF in 
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the distinct model classes. In the second step, these estimated PSFs are em- 
ployed to reconstruct (by deconvolution with the estimated PSFs) a target 
image, and the accuracy of the resulting reconstructions is compared. 

Inference on the PSF as required in the first step has to be conducted 
from dim images, and hence requires low-dimensional modeling. A possible 
approach is to use a parametric model; however, this involves the risk of 
misspecification. As an alternative, one could seek nonparametric estimates 
for the PSF. Due to the dimness of the image, nonparametric smoothing 
algorithms would require a substantial amount of smoothing. Therefore, the 
essential local feature of the PSF, the steep central peak, would be reduced, 
and hence its optical transfer function would be distorted. Thus, as an actual 
estimate of the PSF for the reconstruction process, the Zernike series esti- 
mates or other smoothed estimates should not be used. However, we argue 
that even for a dim image the Zernike estimate with few Zernike moments 
can be used for recovering the global feature of reflectional symmetry. Av- 
eraging along the estimated axes then reduces the noise level in the PSF 
reconstruction, which improves the reconstruction in step 2. 

Specifically, the true PSF in the first step in the simulations consists 
of a bivariate Gaussian density function with full width at half maximum 
[FWHM] of 250 nm along the y-axis and 250/\/2 nm along the x-axis, and 
the bead used to estimate the PSF is assumed to be 50 nm in diameter. 
Moreover, the (true) peak intensity in the image of the bead is ~22, which 
yields a signal-to- noise ratio for the brightest pixels of ~5. 

We use four models in which we estimate the PSF from the available 
(Poisson-distributed) observations. First, we use two parametric models, 
one correctly specified (i.e., the intensities have the shape of a Gaussian 
density with unknown covariance matrix), the other slightly misspecified 
with intensity function proportional to exp(-l/2(g(x, y))^'^^) where q{x, y) = 
{x,y)Ti~^{x,yY' . Both models are estimated by maximum likelihood. Fur- 
ther, we use two nonsmoothed nonparametric estimates. The first simply 
consists of the observed raw data, for the second we average the raw data 
along the two estimated axes of reflectional symmetry, thereby reducing the 
noise level by a factor 2. 

In the second step we aim to recover the target image plotted in Figure 
8 from Poisson-observations with intensities given in (19), that is, the con- 
volution of the target image and the true PSF described above. The target 
image is of size 8.2 pm along the x and y-directions and with 128 x 128 pixels 
along each axis, that is, the resolution of a pixel is ~64 nm. The signal-to- 
noise ratio of the brightest pixels is ~20 (and correspondingly lower for most 
of the image). For the image reconstruction by deconvolution, the distinct 
estimated PSFs are employed in the same algorithm. We use the Expec- 
tation Maximization method [cf. Shepp and Vardi (1982)], also called the 
Richardson-Lucy algorithm [cf. Richardson (1972) and Lucy (1974)], which 
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Fig. 8. Test image used m the simulations of the benefit from using an estimated axis of 
symmetry for the PSF. Left: true image; right: convolved image with Gaussian noise. 



is one of the most commonly used algorithms for deconvolution problems 
with positivity constraint. For each estimate of the PSF we record the small- 
est Li- and -L2-distances attained between any iterate of the Richardson- 
Lucy reconstruction based on the respective PSF and the true target image 
in Figure 8. 

Table 1 shows the mean optimal Li- and L2-distances from 200 simula- 
tions of the imaging process, that is, subsequent execution of steps 1 and 
2. It turns out that while the correctly specified parametric model performs 
best for recovering the target image, symmetrizing the nonparametric es- 
timate greatly improves its performance, even beyond that of the slightly 
misspecified parametric model. 

6. Conclusions. Detection and estimation of symmetry are fundamental 
concepts in many areas of science and technology. In particular, the con- 
cept of symmetry plays an important role in image analysis and pattern 
recognition. 

Symmetry is also relevant in many statistical models. An important and 
well-studied example is the symmetric location model h{x — 9), where h{x) = 



Table 1 

Mean optimal Li- and L'z-distance achieved between reconstructed image and true image, 
based on three different estimates of the PSF 



Distance 


Parametric 


Parametric 


Nonparametric 


Nonparametric 


measure 




(misspecified) 


without symmetry 


with symmetry 


Li (xlO^) 


1.3 


1.9 


2.0 


1.7 


L2 (xW) 


0.6 


1.6 


1.5 


1.2 
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h{—x) is an unknown symmetric density function and G M is the location 
parameter. Such models consisting of a Euclidean parameter as well as a 
nonparametric component are called semiparametric, and efficient, that is, 
asymptotically optimal estimation procedures in such problems are impor- 
tant and difficult issues in statistical inference [Bickel et al. (1993)]. 

In this paper we have discussed how to estimate the angle of the axis 
of reflectional symmetry of an image function, and studied its asymptotic 
properties. This problem is also of a semiparametric form, with the angle (5 
as the target parameter, and the image function (that is reflection symmetric 
with respect to a fixed axis, say, the x-axis) as nonparametric component. 
Although we showed that the parametric rate is achievable for estimating 
the parameter /3, and also obtained an asymptotic normal law, we did not go 
into the problem of semiparametric efficiency and leave this issue for future 
research. 

We have applied our method to calibrating the point-spread function 
(PSF) in confocal microscopy. In particular, we have shown how reflection 
symmetry (but no rotational invariance) may arise in the PSF. Further, we 
demonstrated that estimating the symmetry axes and symmetrizing the im- 
age of the PSF reduces the noise level in nonparametric estimates, and can 
lead to substantial improvement in the performance in subsequent image 
reconstruction algorithms. 

Future research will be directed toward elucidating symmetry information 
and estimation in more complex microscopic setups, in particular, in 3D- 
fluorescence microscopy [e.g., 4PI-microscopy as in Bewersdorf, Schmidt and 
Hell (2006)]. 
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SUPPLEMENTARY MATERIAL 

Estimating bilateral symmetry: Technical details 

(DOI: 10. 1214/10- AOAS343SUPP; .pdf). Here we provide the technical proofs 
for our results in the paper "Improving PSF calibration in confocal micro- 
scopic imaging — estimating and exploiting bilateral symmetry." 
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